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1 What is tight binding? 

"Tight binding" has existed for many years as a convenient and transparent model for the 
description of electronic structure in molecules and solids. It often provides the basis for 
construction of many body theories such as the Hubbard model and the Anderson impurity 
model. Slater and Koster call it the tight binding or "Bloch" method and their historic 
paper provides the systematic procedure for formulating a tight binding model. 1 In their 
paper you will find the famous "Slater-Koster" table that is used to build a tight binding 
hamiltonian. This can also be found reproduced as table 20-1 in Harrison's book and 
this reference is probably the best starting point for learning the tight binding method. 2 
Building a tight binding hamiltonian yourself, by hand, as in Harrison's sections 3-C and 
19-C is certainly the surest way to learn and understand the method. The rewards are very 
great, as I shall attempt to persuade you now. More recent books are the ones by Sutton, 3 
Pettifor 4 and Finnis. 5 In my development here I will most closely follow Finnis. This is 
because whereas in the earlier literature tight binding was regarded as a simple empirical 
scheme for the construction of hamiltonians by placing "atomic-like orbitals" at atomic 
sites and allowing electrons to hop between these through the mediation of "hopping inte- 
grals," it was later realised that the tight binding approximation may be directly deduced as 
a rigorous approximation to the density functional theory. This latter discovery has come 
about largely through the work of Sutton et al. 6 and Foulkes; 7 and it is this approach that 
is adopted in Finnis' book from the outset. 

In the context of atomistic simulation, it can be helpful to distinguish schemes for the 
calculation of interatomic forces as "quantum mechanical," and "non quantum mechani- 
cal." In the former falls clearly the local density approximation (LDA) to density functional 
theory and nowadays it is indeed possible to make molecular dynamics calculations for 
small numbers of atoms and a few picoseconds of time using the LDA. At the other end of 
the scale, classical potentials may be used to simulate millions of atoms for some nanosec- 
onds or more. I like to argue that tight binding is the simplest scheme that is genuinely 
quantum mechanical. Although you will read claims that the "embedded atom method" 
and other schemes are LDA-based, tight binding differs from these in that an explicit cal- 
culation of the electron kinetic energy is attempted either by diagonalising a hamiltonian, 
which is the subject of this lecture; or by finding its Green function matrix elements which 
is the subject of the lecture by Ralf Drautz. 8 The enormous advantage of the latter is that 
calculations scale in the computer linearly with the number of atoms, while diagonalisa- 
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tion is 0(N 3 ). At all events, tight binding is really the cheapest and simplest model that 
can capture the subtleties in bonding that are consequences of the quantum mechanical 
nature of the chemical bond. Some well-known examples of these quantum mechanical 
features are magnetism, negative Cauchy pressures, charge transfer and ionic bonding; and 
of course bond breaking itself which is not allowed by simple molecular mechanics mod- 
els. At the same time tight binding will reveal detailed insight in the nature of the bonds 
and origin of interatomic forces in the system you are studying. 

1.1 The two centre approximation 

In density functional calculations, the hamiltonian is constructed after making a choice of 
functions used to represent the wavefunctions, charge density and potential. If these are 
atom centred, for example gaussians, "fire balls" or Slater type orbitals rather than plane 
waves, then matrix elements of the hamiltonian may become spatial integrals of three such 
functions. An explicit formula taken from the LMTO method is displayed in equation (26) 
in section 3.2 below. This can be the most time consuming part of a bandstructure calcula- 
tion, compared to the subsequent diagonalisation. In the tight binding approximation, we 
side step this procedure and construct the hamiltonian from a parameterised look up table. 
But the underlying theory has the same structure. Each hamiltonian matrix element is con- 
ceived as a integral of three functions, one potential and two orbitals centred at three sites. 
(We have made the Ansatz that the effective potential may be written as a sum of atom cen- 
tred potentials.) If all are on the same site, this is a one centre, or on-site matrix element; if 
the orbitals are on different sites and are "neighbours" while the potential is on one of these 
sites we have a two centre matrix element, or "hopping integral." All other possibilities, 
namely three centre terms and overlap of orbitals on distant sites are neglected. This forms 
a central tenet of the tight binding approximation — the nearest neighbour, two centre ap- 
proximation. The canonical band theory 9 allows us to isolate these terms explicitly and to 
predict under what circumstances these are indeed small (see section 3.2). The two centre 
approximation is more than just a convenient rejection of certain terms; it is implicit in 
the Slater-Koster table and in the calculation of interatomic force that the hamiltonian can 
be written in parameterised two centre form. This allows one to express the dependence 
of hopping integrals upon distance analytically. It is a feature of the quantum mechanical 
method that whereas the hamiltonian comprises short ranged two centre quantities only, 
the solution of the Schrodinger equation using this simple hamiltonian results in a density 
matrix that is possibly long ranged and includes many-atom interactions. Indeed the bond 
order potential exposes this many-atom expansion of the total energy explicitly. 8 

1.2 D(N 3 ) and O(N) implementations 

The obvious way to tackle the tight binding electronic structure problem is the same as 
in density functional theory, namely by direct diagonalisation of the hamiltonian to obtain 
eigenvalues and eigenfunctions in the tight binding representation, section 2.1 below. This 
scales in the computer as the third power of the number of orbitals in the molecule or in 
the unit cell. In the solid state case one employs the Bloch theorem. 10 This means that one 
retains only the number of atoms in the primitive unit cell (rather than an infinite number) 
at the expense of having to diagonalise the hamiltonian at an infinite number of k-points. 
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Luckily there is a well known and sophisticated number of ways to reduce this to a small 
number of points within the irreducible Brillouin zone. 11 ' 12 The Bloch transform of a real 
space matrix Hrl r' l' (in the notation described at equation (3) below) is 

^RLR'L'(k) = ^ %+T)IR'L' e lk T , 
T 

where R and R' run only over atoms in the primitive unit cell, while T are all the trans- 
lation vectors of the lattice. As long as the matrix i?(R + x)L R' v is short ranged this can 
be done easily; for long ranged matrices such as the bare structure constants of (30) below, 
this must be done using the Ewald method. If you like you can define a two centre matrix as 
one for which the Bloch transformation can be reversed (using all A/k points in the whole 
Brillouin zone) 

%+T)lR'L- = -JT ^ -gRLR'L'(k) C~' k T . 
k 

Indeed this is a way to extract a two centre tight binding hamiltonian from an LDA band- 
structure calculation; 13 an alternative approach is described in section 3.2 below. 

In this lecture, I will concentate solely on the method of direct diagonalisation, but 
an alternative and potentially much more powerful approach is to abandon k-space, even 
for a periodic solid, and employ the recursion method to calculate not the eigenvalues and 
eigenfunctions of the hamiltonian H, but its greenian or Green function; formally for a 
complex variable z 

G(z) = (z-H)- 1 . 

Throwing away k-space will lead to a huge computational benefit, namely that the cal- 
culation scales linearly with the number of orbitals, but there is a heavy price to pay — 
interatomic forces converge more slowly than the energy since they require off-diagonal 
greenian matrix elements and the sum rule derived in equation (16) below is not auto- 
matically guaranteed. 14 ' 15 This can play havoc with a molecular dynamics simulation. 
The problem has been solved by the bond order potential which leads to a convergent ex- 
pansion of the tight binding total energy in one-atom, two-atom, three-atom. . . terms — a 
many-atom expansion. This is the subject of the lecture by Ralf Drautz in this workshop. 8 

2 Traditional non self consistent tight binding theory 
2.1 Density operator and density matrix 

The traditional non self consistent tight binding theory, as described, say, by Harrison, 2 
is explained here by following Horsfield et al. 16,17 We use H° to denote the hamiltonian 
to indicate that this is the non self consistent approximation to density functional theory 
as it appears in the Harris-Foulkes functional 5 — the first two lines in equation (36) below. 
(We follow the usual practice of supressing the "hat" on the hamiltonian operator.) Hence, 
H° is the sum of non interacting kinetic energy and the effective potential generated by 
some input, superposition of atom centred, spherical charge densities. 5 The hamiltonian 
possesses a complete set of orthogonal eigenfunctions by virtue of the time independent 
Schrodinger equation, 

H°1p n = £„Vn, 
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which we will write using Dirac's bra-ket notation as 

H a \n) = e n \n) . (1) 
e n are the eigenvalues of the hamiltonian and these are used to construct the band energy, 

Sband, thus 

-Eband = ^ fn £ n - (2) 
n 

Here, /„ are occupation numbers. In an insulator or molecule assuming spin degeneracy 
these are either zero or two depending on whether e„ is greater than or less than the Fermi 
energy. In a metal or molecule having a degenerate highest occupied level these are set 
equal to twice the Fermi function or some other smooth function having a similar shape. 12 
As with any electronic structure scheme, if this is implemented as a bandstructure program 
and hence the hamiltonian is Bloch-transformed into k-space, then the eigenstates are la- 
belled by their band index and wavevector so that in what follows, the index n is to be 
replaced by a composite index nk. (At the same time matrices become complex and you 
may assume that what follows until the end of this subsection applies separately at each 
k-point.) 

Central to the tight binding approximation is the expansion of the eigenstates of H° in 
a linear combination of atomic(-like) orbitals (LCAO). This means that we decorate each 
atomic site, which we denote R to label its position vector with respect to some origin, 
with orbitals having angular momentum L = Im. In this way, I labels the orbitals as s, p 
or d character, while the L label runs as s, x, y, z, xy and so on. These orbitals may be 
written in bra-ket notation as 

|RL) = \i) (3) 

so that we can abbreviate the orbital site and quantum numbers into a single index i or 
j, k, I. In this way we have 

i 

and we use the famous Einstein summation convention, for brevity, whereby a summation 
over the indices i, j, k, I is understood if they appear repeated in a product. (Conversely we 
use n and m to label eigenstates of H in equation (1) and these are not summed implicitly.) 
The expansion coefficients c" are the eigenvectors of H° in the LCAO representation. The 
parameters of the tight binding model are the matrix elements of the hamiltonian in the 
LCAO basis which we write 

We may assume that our chosen orbitals are orthogonal to each other, but to be more 
general there will a matrix of overlap integrals that may also comprise a part of our tight 
binding model. These are 

sa = m ■ 

It then follows from (4) that (summing over j, remember) 

(i\n) = Sijc™ and (n\i) = cjSji (5) 
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in which a "bar" indicates a complex conjugate. The Schrodinger equation (1) becomes a 
linear eigenproblem, 

(H9.-e n S ii )^ = 0. (6) 

In the case of an orthogonal tight binding model, we have Sij = Sij , otherwise we need 
to solve a generalised eigenproblem which is done by a Lowdin transformation. Denoting 
Hfj and Sij in bold by matrices, we insert S~2 S2 after the right parenthesis in (6) and 
multiply left and right by S - ^ : 

= (V^H°S-3 _ £n i) (s'cS - *) = (H - e„l) z, 

which can be solved as an orthogonal eigenproblem, and recover c from z by back- 
substitution using the previously obtained Cholesky decomposition of S. Now we have 
our eigenvectors c™ from which we construct a density matrix, which is central to the 
electronic structure problem. The density matrix provides us with the band energy, local 
"Mulliken" charges, bond charges (in the non orthogonal case) 5 , bond orders, 4 interatomic 
forces, and in the case of time dependent tight binding the bond currents via its imaginary 
part. 18 The density operator p needs to have the following properties. 

Property 1. Idempotency, meaning p 2 = p, 

Property 2. Tr p = N, the number of electrons, 

Property 3. Tr pH° = J2 n /« e « = ^band, the band energy, 

Property 4. Tr pJ^H = Jj-Eband, the Hellmann-Feynman theorem. 

We know from quantum mechanics 19 ' 20 that the one particle density operator is defined as 

P = ^2-fn \n)(n\. 

n 

To find its representation in the LCAO basis, we first define a unit operator, 

i = |i>s« 1 01- (7) 

To show that it is the unit operator, write 

(n\n) = 1 = (n\i) Sr. 1 (j\n) 
= c^SkiS^ Sjicf 

— n n 

— Cj OijCj 

(after using (5) and swapping indices) which is consistent with (4). More generally we 
have 

(n\m) = S nm = c?S ijC f. (8) 

Now using our unit vector, we write the density operator in our, possibly non orthogonal, 
LCAO basis, 

P = £/n l«> (n| = ]T/„i|n) (n|l 

n n 

= J2fn \i)c^{j\. (9) 

n 
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A matrix element of the density operator is 

P« = X> (k\i)c?c-(j\l) 

n 

= Y,fnS ki C^S jl (10) 

n 

and in an orthogonal basis this reduces to the familiar density matrix 

n 

If you are familiar with general relativity or non cubic crystallography then you may wish 
to view the matrix SV, as the metric tensor that "raises" and "lowers" indices of covari- 
ant and contravariant vectors. 6 ' 15 ' 21 Finnis 5 makes this point by distinguishing between 
"expansion coefficients" and "matrix elements" of the density operator. In this way the 
expansion coefficients of the density operator in the LCAO basis are J2 n fn c ?Cj> while to 
obtain density matrix elements their indices are "raised" by elements of the metric tensor 
as in (10); in the orthogonal case (SV, = Sij) this distinction vanishes. 

Now we can demonstrate that p has the properties 1-4 above. The following is really 
included here for completeness as the student may not find it elsewhere in the literature. 
However, on a first reading you may skip to section 2.3 after looking at equations (11), 
(12), (13), (16) and (17). 

Property 1. Idempotency follows immediately from (9). 
Property 2. Tr p = N. We must take the trace in the eigenstate basis, hence 
Tr ' 5 = EE/« Hi)c n c-(j\m) 

m n 
m n 

^ ^ ^ ^ fn ^mn^nm ^ ^ fn 
m n n 

After the second line we have used (8). We can make partial, "Mulliken" charges qi 
which amount to the occupancy of orbital i, 

i n 

using (8). Because of its importance in tight binding, we will write the Mulliken 
charge associated with orbital i explicitly, 

n j 

in which the sum over i implied by the summation convention is, in this instance, 
supressed. This is a weighted decomposition of the norm. Note that in this and the 
following you can easily extract the simpler expressions for the more usual orthogonal 



6 



tight binding by replacing Sij with the Kronecker 5ij in the implicit sums, in which 
case 

n 

It is worthwhile to note that in an orthogonal tight binding model the total charge can 
be decomposed into individual atom centred contributions; on the other hand non or- 
thogonality introduces bond charge 4 so that as seen in (1 1) there is a summation over 
both atom centred and bond charges. You may prefer the latter picture: we all know 
that in a density functional picture the covalent bond arises from the accumulation of 
charge in between the atoms; in an orthogonal tight binding model one might ask how 
is this accumulation described? The answer is that it is captured in the bond order. 

Property 3. Tr pH° = J2 n fn £ n = E h&nd . 

Tr pH° = £ £ /„ (m\i) c?c? 01 H° \m) 

m n 
m n 

^ ^ ^ ^ fn ^mn^nm^m ^ ^ fn -^band 
m n n 

using (1). One may wish to construct partial band energies, E i7 in an equivalent way 

as 

#band = ^2 Ei = fn 'BijCj ■ 

i n 

The corresponding decomposition of the bond energy (18) in section 2.3 is the starting 
point of the many-atom expansion in the bond order potential. 8 

Property 4. The Hellmann-Feynman theorem tells us that 

because solution of the eigenproblem (6), through the Rayleigh-Ritz procedure leads 
us to a density matrix that is variational with respect to any parameter A which may 
be, for example, a component of the position vector of an atom R. Hence to calculate 
the interatomic force we need to find 

Tr P§x H ° = E E /« Hi> 01 Y X H* \m) 

d 

Tx 1 



Now our tight binding model furnishes us with hopping integrals, i??-, and by em- 
ploying a suitable scaling law, for example equation (23) below, the two centre ap- 
proximation and the Slater-Koster table we will know how these depend on bond 
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lengths and angles; so while we don't actually know (i\ -§^H° \j), the derivatives that 



we do know are 



'<9A J 



<9A 



So 



3A a A y <9A 

n L 

Now, to deal with the unknown last two terms, using (4) 

E^s" 



= E/« 



<9 5 
c?(— i|n)e„ + £„(n|— j)cj 



1 6>A- 



d 



since 



Finally we arrive at 



d 



Tr %tf° = E« c ? 



dX » n dX ij 



(12) 



2.2 Density of states and bond order 

The density of states is central to electronic structure theory and is defined to be 22 

n(e) = ^ *(£-£„). (13) 

n 

We can define a partial or /oca/ density of states, nj(e), which is the density of states 
projected onto the orbital i. We write 

n(e) = ^(n\S(e-H°)\n) 

n 

= EEW^(J'H (m|,5( e -ff )|n> 

n m 

-EE tjSjiSjSjkd? (m| % - |n) 

n m 

= Y,^S jk cl (n\5(e-H°)\n). 

n 

The first line follows from the Schrodinger equation (1) and in the second line we have 
inserted our unit operator (7) and a further unit operator, J2 m \ m ) ( m \- The fourth 
line follows because of the orthogonality of the eigenvectors, \n) which means we have 
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(m\S(e — H°) \n) = (n\ 5{e — H°) \n) 8 mn . Remember that in the fourth line j and k are 
dummy orbital indices to be summed over. We can replace these with i and j for neatness 
and this leads to 

n(e) =J2<%S ijC ]6(e - e n ) = ^^(e). (14) 

n i 

Writing the summation over j explicitly we see that the local density of states is 

^( £ ) = EE^-^(e-e„), d5) 



with no summation over i, and that this is a weighted density of states; the weight in an 
orthogonal basis is simply |c™ | 2 — compare this with the Mulliken decomposition (11). 

An example is shown in figure 1. This is a self consistent magnetic tight binding 
calculation of the electronic structure of a Cr impurity in Fe, modelled as a dilute, ordered 
FeisCr alloy. 23 Very briefly magnetic tight binding is achieved by including a spin index, 
|i) = \HLa), (now the occupation numbers vary between zero and one, not two) and 
adding an exchange potential to the self consistent hamiltonian to allow these to split. In 
addition to the Hubbard-?/ (see section 4) one includes a "Stoner I" parameter. We cannot 
go into details here, but it's gratifying that the simple tight binding model quantitatively 
reproduces the LSDA result, even to the extent of predicting the "virtual bound state" on 
the Cr impurity. 24 ' 25 

The density of states can be used to find the band energy, since by the properties of the 
Dirac delta function, 

E f" / S (S - £n)ede = /« £ ™ = ^band- 

n n 

If we allow the occupation numbers to be represented by the spin degenerate Fermi-Dirac 
distribution, 2/(e), then we find, using (13) and our property 3, above, 

£band = 2 / f(e) e n(e) de = Tr pH° (16) 



which is an important identity in tight binding theory and one which bears heavily on the 
convergence of the many atom expansion in the bond order potential. 26 

Finally in this section we should mention that the bond order which is central to the 
bond order potential 8 is obtained directly from the density matrix elements. We define 

@ij = ^ (Pij+Pji) 

as the partial order of the bond as contributed by orbitals i and j, it being understood that 
these are on different atomic sites. The bond order between sites R and R' is obtained by 
summing the partial bond order over all the orbitals on each atom in question, 

e RI r = ]Te RLR ^. (17) 

LL' 
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Figure 1. Example of a local density of states. 23 This is an ordered alloy, FeisCr on a body centred cubic (bcc) 
lattice. On the left is the local spin density functional result and on the right a simple, non orthogonal magnetic 
tight binding approximation. As is conventional, the spin up and down densities are shown as upright and upside 
down functions respectively. The Fe atom shown is the one closest to the Cr impurity and the density is projected 
onto the c(-manifolds. Apart from the accurate description provided by the tight binding model, the most striking 
feature is the virtual bound state 24,25 seen as sharp peak in the local Cr density of states. It's notable that the 
occupied, spin up state has t^g symmetry while its unoccupied partner belongs largely to the e g manifold. 



2.3 The tight binding bond model 

Just as in density functional theory, the sum of occupied eigenvalues of the one electron 
hamiltonian is not the total energy. In the traditional tight binding approximation, begin- 
ning probably with the papers of Jim Chadi, 27 one writes simply 

Etot = -Eband + -Epair 

for the total energy in the band model and E pa i r is a pairwise repulsive energy whose func- 
tional form and parameters constitute ingredients of the tight binding model; it is intended 
to represent the double counting and ion-ion contributions to the density functional total 
energy. 27 "Double counting" is a term given to the electron-electron interaction energy 
in density functional theory. Because the theory is cast into a one electron form through 
the Kohn-Sham equations, the band energy, by summing over the eigenvalues, counts the 
electron-electron interaction twice. The interaction between, say, electrons in occupied 
states 1 and 2 is counted first when eigenvalue 1 is added in and again when eigenvalue 2 is 
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added. One cannot simply divide by two because Eb an d also contains kinetic and electron- 
ion energies which are not double counted. Hence one recalculates the electron-electron 
interaction energy and subtracts it, calling this the "double counting" correction. 

Pursuing an argument that goes back as far as the sixties, 28 ' 29 Pettifor 30 formulates the 
total energy in terms of the bond energy, Eb on d, rather than the band energy. The tight 
binding bond model 6 (TBBM) is the starting point for both self consistent tight binding 
which is described below in section 4 and for the modern bond order potentials. 8 Therefore 
we will pursue only the bond model further here. The essential point is that one arrives at 
the covalent bond energy*' 6 by removing the diagonal elements of Eb an d = Tr pH °. We 
recall that orbital indices i and j are a composite of site labels and quantum numbers, and 
write 

2 X! 2 PRLR>L> #R/L'RL- (18) 

RI R'L' 

Here all terms are excluded from the double sum if orbitals i and j are on the same site 
R. Note how by dividing and multiplying by two we can expose this as a sum of bond 
energies which is then divided by two to prevent each bond being double counted in the 
same way as a pair potential is usually written. 

In the TBBM, the remaining diagonal terms in Eb an d are grouped with the correspond- 
ing quantities in the free atom. In the non self consistent tight binding approximation, the 
on-site matrix elements of H° are simply the free atom orbital energies (eigenvalues of the 
atomic hamiltonian) 

and in addition to the hopping integrals, these are parameters of the tight binding model, 
e s , s p and Ed- Furthermore, we assume certain orbital occupancies in the free atom, say, 
Nm, whereas after diagonalisation of the tight binding hamiltonian one finds these orbitals 
have occupancy given by the diagonal matrix elements of the density matrix. Hence there 
is a change in energy in going from the free atom limit to the condensed matter which is 

E prom = (PRLRL ^RLRi - ^R£ £Rl) 
RL 

= ^2 (prlrl - Nw) £ Ri 

RL 

= J2 A lRLe Re . (19) 

RL 

We have assumed for now that on-site elements of H° are strictly diagonal and we recog- 
nise the first term in the first line as the difference between Eb an d and -Ebond- E prom is 
called the promotion energy since it is the energy cost in promoting electrons that is very 
familiar, say, in the s-p promotion in covalent semiconductors in "preparing" the atoms in 
readiness to form the sp 3 hybrids in the diamond structure or the sp 2 hybrids in graphite. 
Thus in the tight binding bond model the binding energy is written as the total energy take 
away the energy of the free atoms, 

Eq = Ebond + E prom + Ep a i r . (20) 



Ebond - j X! 2 P*i B % = 

ij 
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The interatomic force is minus the gradient of the pairwise i? pa i r which is trivial, minus 
Tr pVi? which can be computed using equation (12) assuming that on-site hamiltonian 
matrix elements remain constant; this is the fully non self consistent tight binding ap- 
proximation. And in fact at this level of approximation the band and bond models are 
indistinguishable. The first order variation of Eb with respect to atom cartesian coordinate 
R a is 

J^E B =J2 2 PRiR ^^-^, i , RL + ^-i? prom + ^-i? pair . (21) 

R.VR, 

This is written for the orthogonal case, since this approximation forms a tenet of the 
TBBM. However, it's easy enough to add in the term from (12) containing the overlap 
and of course the diagonal elements Su are constant and do not contribute to the force. 
Note that the half in front of (18) has vanished — in the calculation of the force one sums 
over all bonds emanating from the atom at R, not just half of them! 

Now comes a rather subtle point. Unlike the band model, the bond model is properly 
consistent with the force theorem. 31 This states that there is no contribution to the force 
from self consistent redistribution of charge as a result of the virtual displacement of an 
atom. If a self consistent electronic system is perturbed to first order then that change in 
the bandstructure energy due to electron-electron interaction is exactly cancelled by the 
change in the double counting. This remarkable result means that by making a first order 
perturbation one cannot distinguish between an interacting and a non interacting electron 
system. 32 Indeed to calculate the interatomic force it is sufficient to find the change in 
band energy while making the perturbation — in this case the virtual displacement of an 
atom — in the frozen potential of the unperturbed system. In the band model there will be 
a first order change in the band energy upon moving an atom which ought to be cancelled 
by an appropriate change in the double counting, but is not because this is represented by 
the pair potential. Now we can discuss dE pmm /dR a . In the band model there is no con- 
tribution to the force from -E prom (19); because of the variational principle Crl^rl = 0, 
and qRLSeRL = because the erl are constants. However the Mulliken charge transfers 
are not necessarily zero and the force theorem does require any electrostatic contributions 
due to charge transfer to be included in the interatomic force; 33 ' 34 neglect of these leads 
to the inconsistency of the band model. In the TBBM the most limited self consistency is 
imposed, namely the Ansatz of local charge neutrality so that electrostatic charge transfer 
terms vanish. This requires that for each site the total Mulliken charge difference between 
free atoms and condensed phase summed over all orbitals is zero. This is achieved iter- 
atively by adjusting the on-site orbital energies. Here is the simplest example of a self 
consistent tight binding theory. It only affects the diagonal, on-site hamiltonian matrix 
elements and hence only -Eprom is changed. Suppose we now write the hamiltonian as 

H = H° +H' (22) 

where H' has only diagonal elements which we may call lS.e-R.L- Then 

RL 

= ^2 Ag RL (e R£ + Ae Re ) . 

RL 
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In a sense this isn't really "promotion energy" anymore because we have applied the on-site 
energy shift to the free atoms also, but it is consistent with the formulation of the TBBM. 6 
There will now be a contribution to the force on atom R from the new term Aq^ Aeg. 
If the self consistency is achieved in such a way that all orbital energies are shifted by the 
same amount at each site, then this contribution vanishes because Asi is independent of L, 
moves to the front of the summation sign and ^2 L Aql = by the local charge neutrality 
condition. Further and complete discussion of the TBBM can be found in the original 
paper 6 and in Finnis' book. 5 

3 How to find parameters 

Now we turn to the question that is probably the most controversial. Many people dislike 
the tight binding approximation because whereas on the one hand we claim it to be close 
to the ab initio local density approximation solution, on the other we are reduced to finding 
parameters empirically just as if this were another classical potential. My own view is 
that if the tight binding approximation contains enough of the physics of the system we 
are studying then any reasonably chosen set of parameters will provide us with a useful 
model. From this point of view we would also demand that only a very small number 
of parameters is actually employed in the model. Furthermore it should be possible to 
choose these by intelligent guesswork and refinement starting from some well established 
set of rules; for example Harrison's solid state table, 2 or the prescription of Spanjaard and 
Desjonqueres for the transition metals. 35 For example, the latter prescription has furnished 
us with useful tight binding models 23 ' 36 for Mo, Re, Nb and Fe each with some five to ten 
adjustable parameters. Alternatively a 53-parameter model for Mo was produced by very 
careful fitting to a huge database of properties. 37 There doesn't appear to exist a particular 
advantage of one approach over the other and both types of model have turned out to be 
predictive of electronic and structural properties of the transition metals. 

We need to distinguish between hamiltonian parameters — on-site orbital energies ejn 
and hopping integrals H RLR , L , — and the parameters of the pair potential. Additional com- 
plications arise as described later in section 3.3 in the case of environmentally dependent 
parameters. 37 

I wish to illustrate the problem by reference to some examples from the literature. 
3.1 Parameters by "adjustment" — example of Zr02 

The tight binding model for zirconia, 38 ZrC>2, was designed to provide a description of the 
structural properties of this industrially important ceramic material. ZrC>2 suffers a num- 
ber of structural phase transitions as a function of temperature. This is exploited in an 
extraordinary phenomenon called transformation toughening. 39 Its low temperature phase 
is monoclinic, at intermediate temperatures it is tetragonal and the high temperature mod- 
ification is cubic. An open question was whether the tetragonal to cubic transition is of 
first or second order thermodynamically, order-disorder or displacive. Additionally, it is 
known that the cubic structure is stabilised at low temperature by doping with aliovalent 
cations (Y, Ca, Mg etc) while the mechanism for this was unknown. The tight binding 
model turned out to be capable of addressing both these issues and the order of the transi- 
tion was discovered 40 as well as the mechanism of stabilisation of the cubic phase. 41 The 
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Figure 2. Bond integrals, after Majewski and Vogl. 42 This shows the well known atomic orbitals of the various 
s,pord types joined along a bond. Radial symmetry along the bond is assumed leading to the designation of the 
bond as a, n or 5. To construct a tight binding hamiltonian requires these fundamental bond integrals assembled 
through the Slater-Koster table using the direction cosines of the bond in a global cartesian system (these bond 
integrals are given with respect to a z-axis directed along the bond). This is illustrated in figure 6.5 in ref [3]. 



strategy of finding tight binding parameters was quite simple. Since the eigenvalues of 
the hamiltonian describe the energy bands it is sensible to adjust the on-site energies and 
hopping integrals to the LDA bandstructure, and then find a simple pair potential whose 
parameters are chosen to obtain, say, the equilibrium lattice constant and bulk modulus. In 
this case the smallest number of adjustable parameters was chosen to replicate the cubic 
phase in the hope that the model will then predict the ordering in energy of the competing 
phases. The steps are these. 

1 . Choose a minimal tight binding basis set. In this case e?-orbitals were placed on the 
Zr atoms and s and p on the oxygen. We should mention that being an ionic crystal 
the TBBM is inadequate and this is in fact a self consistent tight binding model using 
polarisable ions. This is explained later in section 4. The hopping matrix elements are 
linear combinations of the fundamental bond integrals that are illustrated in figure 2. 
The particular linear combination depends on the bond angle geometry and is encap- 
sulated in the Slater-Koster table. 1 This is illustrated in figure 6.5 in ref [3]. We only 



14 



Fluorite (LDA) Fluorite (TB) 
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Rutile (LDA) Rutile (TB) 




Figure 3. Energy bands of Zr02 using both LDA and a tight binding model in both fluorite and rutile crystal 
modifications. The model parameters were adjusted to the fluorite bands and the rutile bands are therefore a 
prediction. We should note that a number of features such as the splitting in the c(-manifold into *2 and e s sub- 
bands and the crystal field widening of the p-derived ligand band in rutile are consequences of using the self 
consistent polarisable ion model, and this will be described later in section 4. But we can note in anticipation 
that it is the new A parameters that permit the ordering (t'2 > e g ) in the cubic crystal field and vice versa in the 
octahedral field to be reproduced automatically. 



need to find the relevant fundamental bond integrals between neighbouring atoms. 
Zr-O first neighbour bonds require us to know sda, pda and pdir and we choose also 
to include second neighbour O-O bonds to be made by ssa, spa, ppa and ppir bond 
integrals. We have to choose both their value and the way in which they depend on 
bond length. There is a "canonical band theory," that is really appropriate for met- 
als 9 ' 43 ' 44 but which faux de mieux we can apply more generally. This provides us with 
guidance on how the bond integrals decay with distance and also with certain ratios, 
namely ppa:ppir and dd<r:ddir:dd5, see equation (30) below. The required hopping 
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Figure 4. Total energy versus volume in four competing crystal structures of ZrC>2 . 38 At each volume, the energy 
is minimised simultaneously with respect to all the remaining degrees of freedom, (a) LDA calculations of the 
absolute binding energy (energy with respect to spin polarised free atoms); (b) tight binding results referred to 
the equilibrium energy of the monoclinic phase, (c) and (d) show the axial ratio q and distortion parameter 5 in 
the tetragonal modification as a function of volume. 



integrals are initially taken from Harrison's solid state table and adjusted visually to 
obtain agreement with the shapes, and especially the widths of the LDA bands. One 
can also adjust to either the LDA or to experimental band gaps. Also the scaling of the 
bond integrals can be adjusted to the volume dependence of the LDA bandwidths. a 
The result is shown in figure 3. 

We should give more detail of how the bond integrals depend on bond length, r. A 
very useful function is that of Goodwin, Skinner and Pettifor 46 (GSP) 



{(.I'm) = Vb I - J 'Mi \ n 



(23) 



a It is very useful to have a computer program that can calculate energy bands, density of states, total energy using 
both LDA in some form and in the tight binding approximation, preferrably all using the same input file. Luckily 
such a program exists. 45 Students may contact the author if they wish to learn how to use this. 
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Most important are the prefactor Vq which is the value at the equilibrium bond length, 
d, and the exponent n which determines the slope of the function at equilibrium, since 
when r — d the argument of the exponential vanishes. The role of n c and r c is to give 
a rapid decay to the function at around r = r c . 

2. A pair potential needs to be chosen. The GSP function can be used but in the Zr0 2 
model a very simple Born-Mayer form was used between first neighbour Zr-O bonds 
only. The Born-Mayer function ip(r) = Ae~ br has only two parameters which were 
fitted to the lattice constant and bulk modulus of cubic Zr02- 

Figure 4 shows energy volume curves for the competing crystal structures comparing 
the tight binding model to LDA. Also shown are the order parameters that describe the 
tetragonal to cubic phase transition as functions of volume. 

It is rather clear that the tight binding model for Zr02 gives a really excellent set of 
predictions, having been fitted (or adjusted, rather) only to the cubic structure. In particu- 
lar the rutile structure is found to be much higher in energy than its competitors — a feature 
that cannot be reproduced in purely classical models. The vanishing of the order parame- 
ters with pressure is well reproduced qualitatively. This and the example shown in figure 1 
where simple models, rather insensitive to the choice of parameters, reveal useful and pre- 
dictive physics gives one confidence the tight binding approximation is indeed a valuable 
and reliable theory. 

3.2 Parameters taken from first principles tight binding — example of Mo 

Students who are not particularly interested in the details of an LMTO calculation, may 
skip this section after looking at figure 5 and subsequent comments. However section 3.3 
is important. It makes sense to obtain the hamiltonian matrix elements from ab initio band- 
structures. Probably the most transparent LDA bandstructure theory is the one provided 
by the linear muffin-tin orbitals (LMTO) method. In the atomic spheres approximation 
(ASA) the entire bandstructure problem is reduced to knowing four "potential parameters" 
in each Ri site and angular momentum channel. Moreover these parameters have a clear 
interpretation in terms of the bandstructure. C is the centre of the band; A is the bandwidth 
parameter; 7 is a distortion parameter describing the deviation from canonical bands and 
finally p is a small parameter allowing the eigenvalues to be correct up to third order in 
their deviation from some chosen energy called s v . An LMTO is a composite orbital-like 
basis function. A sphere is enscribed about each atom with radius such that the sum of all 
sphere volumes equals the total volume; in a simple monatomic crystal this is the Wigner- 
Seitz radius. Within the sphere the radial Schrodinger equation is solved at the energy 
e v in the current potential and this solution and its energy derivative are matched to solid 
Hankel and Bessel functions between the spheres. This matching condition is enough to 
provide the potential parameters which are functions of the logarithmic derivatives of the 
radial Schrodinger equation solutions </>i(r) = <j>i(r) Yl{y). Each LMTO envelope may 
be expanded about a given atomic site using the property that a Hankel function at one 
site may be written as a linear combination of Bessel functions at some other site. This 
property means that all the Hankel functions in the solid can be expressed as a "one cen- 
tre" expansion about any one atomic sphere. The expansion coefficients are called "k = 
structure constants" and they transform under rotations according to the Slater-Koster table 
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and hence may be identified as Wm hopping integrals. 47 48 However conventional struc- 
ture constants are very long ranged. To make contact with tight binding theory Andersen 
and Jepsen showed that one can make similarity transformations between sets of solid state 
LMTO's; 49 each basis set being equivalent to another since they give identical bandstruc- 
tures. In particular Andersen demonstrated that one can defined a "most localised" and an 
"orthogonal" set of LMTOs. The transformation works like this. In the ASA an LMTO 
at site R is made up of a linear combination of a radial solution <p(r — R) (the "head") 
and energy derivative functions <p(r — R') (d</>/de evaluated at e v ) at all other sites (the 
"tails"). These are assembled into a one centre expansion using the structure constants. So 
an LMTO looks like this, 

XR L (r-R) =cj> RL (r-R)+ <fi R , L ,(r - R') h WL , RL . 

R'L' 

By a choice of normalisation, one can choose the <j>(r — R') to be those that are orthogonal 
to the radial solutions in each sphere. This particular set of energy derivative functions is 
given a superscript 7 and one is said to be using the "7-representation." More generally 
one can vary the normalisation by mixing in some radial solutions with the <p(r — R') to 
make up the tails of the LMTO. To do this we write 

fet (r — R) = <& L (r - R) + <t> RL (r - R) o RL , (24) 

so that in the 7-representation, the potential parameter ori is zero. It's called o for overlap 
but has units of energy -1 . To construct the overlap matrix in the ASA one has to expand 
out (x|x); and similarly (x \— V 2 + V e g\ x) for the hamiltonian. If we write that /ir/z/ rl 
is an element of a matrix h and orl and p R L are elements of diagonal potential parameter 
matrices, o and p, then Andersen finds for the overlap matrix 48 

S = 1 + oh + ho + hph. (25) 

As we mentioned p is a small potential parameter. So in the 7-representation o — and to 
second order the overlap is unity and we have an orthogonal basis. The hamiltonian matrix 
turns out to be 48 

H = e v + h + hoe„ + £„oh + h (o + pe u ) h. (26) 

Again, in the 7-representation, neglecting third order terms the hamiltonian is just H = 
e v + h. So if one calculates structure constants and self consistent potential parameters 
using an LMTO code then one can build an orthogonal tight binding model by explicitly 
building H in the 7-representation. By construction, to second order it will reproduce the 
LDA energy bands. 

Unfortunately there is no guarantee that this hamiltonian is short ranged. Andersen 
made a particular choice of the potential parameter orl by defining "screening constants" 
«rl in this way: ref [9], eq (91), 

= Crl — £i/,RL — • (27) 

orl Jul - a R L 

They are called screening constants because the effect of adding radial solutions to the 7 
in (24) is to match the Schrodinger equation solutions in the sphere to Hankel functions 
Kl{t — R) that have been screened by additional Hankel functions at surrounding atomic 
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sites. There is an electrostatic analogy. The solid Hankel function represents the elec- 
trostatic potential due to a 2 l multipole. This can be screened by surrounding the sphere 
with further grounded metal spheres, whose contribution to the potential is then provided 
by these further Hankel functions at the surrounding spheres. If one chooses the screen- 
ing constants orl to equal the band distortion parameters 7rl then one arrives at the 
7-representation since we get orx = in (27). All other representations are specified by 
choices of screening constants. The choice qrl = corresponds to the so called "first 
generation" LMTO which employs the standard k = KKR structure constants^ 

- -^(- i ) \ 2e % K2 f: m ^"( r - r ') (2«) 

where 

K L {v)=r- l - 1 Y L {v), 

is the solid Hankel function, 

CWl = JjdriY L „Y L ,Y L (29) 

are Gaunt coefficients and Yj, are real spherical harmonics (see Appendix). The whole 
Slater-Koster table is encapsulated in this formula; the Gaunt coefficents provide selection 
rules that pick out certain powers of r and angular dependencies. By pointing a bond 
along the z-axis one can see how the canonical scaling and ratios come about since these 
structure constants are simply, 48 

B SS a = —2/d 

B spa = 2V3/d 2 
B pp{a ^=6{2,-l}/d 3 

B sda = -2V5/d 3 

B pd{t7 , n} = 6\/5{-\/3, l}/d 4 

B dd{<w} = 10{-6,4,-l}/d 5 (30) 

in which d is a dimensionless bond length r/s, where s is conventionally chosen to be the 
Wigner-Seitz radius of the lattice. These can be compared with the cartoons in figure 2 
in which the overlapping of two positive lobes leads to a negative bond integral and vice 
versa. This is because the orbitals are interacting with an attractive, negative, potential 
(section 1.1). Note how the factor (— I) 1 in (28) neatly takes care of the cases like psa — 
—spa. You have to be careful of these if you program the Slater-Koster table by hand. 5 

Transformations from the "first generation" to "second generation" LMTO basis sets 
are quite easily done. Having chosen screening constants one transforms the structure 
constants thus, c 

B a =B + BaB Q (31) 



''Andersen uses the symbol S for structure constants but we've already used it for the overlap, which is standard 
tight binding usage. Here we use B for Andersen's which differ by a prefactor 2/[(2£ — 1)\\(2£' — 1)!!] and a 
minus sign from the KKR structure constants. 50 

c Clearly B r i l i RL has two centre form, section 1.1, as it depends only on the connecting vector R — R' (28). 
It's less obvious that B Q is a two centre matrix because of the three centre terms introduced by the second term 
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which is a Dyson equation, and a is a diagonal matrix. Then one transforms the potential 
parameters by defining a vector (we supress the RL subscripts) 

f = l + «7- e „)£ 

after which (ref [9], p. 88) 

c = e v + {C-e v )£ ; d = £ 2 A 

where c and cZ are the band parameters C and A in the new representation. The overlap 
parameter o is transformed according to (27). 

Andersen and Jepsen 49 determined empirically a set of screening constants, namely^ 

a s = 0.3485 a p = 0.05304 a d = 0.010714, (32) 

which lead to the "most localised" or most tight binding LMTO basis. Now one can con- 
struct hamiltonians and overlaps according to (26) and (25) by noting that the first order 
hamiltonian is constructed from potential parameters and structure constants 9 ' 48 

hnLR'L' = (crl — £i/,Rl) SrLR'L' + \/ ^RL -BrlR'L' V ^R'L'- 

Now we want our tight binding hamiltonian to have two centre form and it is easy to 
identify which are the three centre terms in the LMTO hamiltonian and overlap matrices — 
they are contained in the terms bilinear in h, the last terms in (26) and (25). These terms 
(as do the linear terms) also contain two and one centre terms, of course, arising from the 
diagonal terms of h. We can dispose of three centre terms in two ways. 

1. We can work to first order, in which case, in both a- and 7-representations 

H (1) =e„ + h (33) 

and since oh terms are of second order, both these are orthogonal models with overlap 
being unity. 

2. We can work to second order by retaining oh terms but neglecting the small potential 
parameter p 1 in the 7-representation. In this representation (o = 0) this is no differ- 
ent from the first order hamiltonian, and the overlap is unity. In the a-representation 
this introduces some additional two centre contributions to the matrix elements of the 
hamiltonian and overlap, and we are careful to extract one and two centre contribu- 
tions from the last term in (26). 

All this is illustrated in figure 5 for the bcc transition metal Mo. The screening constants 
from (32) are used. Here are some noteworthy points. 

1 . Clearly the two representations deliver different sets of hopping integrals. You cannot 
expect density functional theory to furnish you with THE tight binding model. On 
the other hand they show a proper decay with increasing bond length. The decay is 

in (31). Nonetheless because the transformation is done in real space it is also a two centre matrix by virtue again 
of its dependence only upon R — R'. On the other hand it possesses additional "environmental dependence," see 
section 3.3. 

^An alternative is to define a Re = (2£ + l)(r R ^ / s) 2l+l by choosing site and <?-dependent "hard core radii" 
rji(. sl This is consistent with "third generation LMTO." 52 
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Figure 5. Hopping integrals ti'm in the body centred cubic transition metal Mo, calculated using LMTO theory. 
The three integrals dda, ddn, and dd8 are found by rotating the z axis to first and then to second neighbour 
bonds and doing this at three different atomic volumes; 53 hence for each integral six values of tt'm are shown 
as a function of bond length. Three model LMTO hamiltonians are used. The crosses refer to the two centre 
7-representation; the circles to the first order a-representation and the pluses to the second order, two centre 
H" . In the lower panel are shown the diagonal matrix elements and their, rather strong, volume dependence. 



more rapid in the tight binding, a-representation as expected, furthermore the first 
order tight binding representation is strictly orthogonal; not shown in figure 5 are the 
overlap matrix elements in the second order tight binding representation, but indeed 
these are very small — no greater than 0.025 in magnitude. Note that the tight binding 
bond integrals respect the signs and roughly the canonical ratios of the bare structure 
constants (30) while in the 7-representation ddS and the second neighbour ddir have 
the "wrong" signs. Furthermore we would find that while the tight binding bond 
integrals shown reproduce the LDA bands using just first and second neighbour matrix 
elements, this is not the case for the 7-representation. Note that the first and second 
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order tight binding matrix elements are essentially the same; the additional second 
order terms may be safely neglected and the first order orthogonal hamiltonian (33) is 
clearly the proper one to use for this case. 

2. If you have the patience, then you can do this exercise by hand in the case of the first 
order hamiltonian. 48 However the scheme has been automated and is implemented in 
our LMTO suite. 45 

3. Unfortunately the on-site energies shown in the lower panel of figure 5 are far from 
independent of volume. This is a remaining unsolved question for the construction 
of tight binding models in which the on-site energies are invariably constant (except 
of course for the adjustments in self consistent models to account for electrostatic 
shifts due to charge transfer, see (43) below). Andersen 48 points out that the Dyson 
equation (31) provides guidance on how to account for this volume dependence in 
terms of the local neighbour environment. Whereas on-site matrix elements of the 
bare structure constants, B are zero, we have from (31) 

B RLRL= ^2 ^2 B RLR"L" an"i" Br"I"RL 
R"#R L" 

and the on-site matrix element of (33) is 51 

£rl = crl + d RL B^ LRL . 

However the band centre parameter c and bandwidth parameter d are also strongly 
volume dependent. 9 ' 44 An important contrast with the ASA is that in tight binding, 
the on-site parameters are constant — the scaling law has to take care of both the bond 
length dependence at constant volume and the volume dependence itself. 54 



3.3 Environmentally dependent tight binding matrix elements 

Possibly the most striking feature displayed in figure 5 is a discontinuity, most notably in 
the ddir and ddS bond integrals, between first and second neighbours. This is of particular 
importance to structures like bcc which have first and second neighbours rather similar in 
bond length. It means that one cannot find a simple scaling law, such as the GSP (23) 
that can connect all the points in the graph. This effect was noticed in the case of the ssa 
bond integral in Mo by Haas et al? 1 and they proposed a very significant development in 
tight binding theory, namely the use of environmentally dependent bond integrals. 55 The 
discontinuities in the dd bond integrals were noticed by Nguyen-Manh et al. 53 who offered 
the physical explanation in terms of "screening." The basic idea is that the bond between 
two atoms is weakened by the presence of a third atom. Therefore the scaling of a bond 
integral, say by the GSP function (23) is modified by multiplying it by (1 — Sw m ) where 
the "screening function," Su'm, is the hyperbolic tangent of a function 37 

IR-R/'I + IR' -R"i\""'- 



exp 



R," 



IR-R' 



(34) 



in which A, A and 77 are parameters to be fitted. This complicated expression can be simply 
explained. 37 ' 53 As a third atom, R" approaches the R R' bond the term in parenthe- 
ses becomes small, and approaches one in the limit that atom R" sits inside the R R' 
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bond. This increases the value of the exponential and the tanh function smoothly reduces 
the R — R' bond integral. Whereas Tang et al. 55 introduced this function empirically, 
Nguyen-Manh et al. 53 were able to derive its form using the theory of bond order poten- 
tials, and explain why dda is not strongly screened while ddn and dd8 are. Modern tight 
binding models 51 ' 56 ' 57 for transition metals are now fitted to curves such as those in figure 5 
using (34). Indeed in these new schemes a repulsive energy is also fitted to an evironmen- 
tally dependent function similar to (34). This is intended to make a better description of 
the valence-core overlap 44 58 between atoms which is short ranged but not pairwise and is 
otherwise not properly captured in the tight binding bond model. So nowadays one finds 
instead of (20) 

Eb = -^bond + E pIom + E cnv + -Ep a i r (35) 

in the TBBM, and -E onv is the new environmentally dependent repulsive energy; it being 
understood that -Ebond m ^y be constructed using environmentally dependent hopping in- 
tegrals too. Eprom is sometimes omitted, 56 ' 57 in the instance that only one orbital angular 
momentum is included in the hamiltonian, for example if one employs a d-band model for 
transition metals. 



4 Self consistent tight binding 



We described a tight binding model for ZrC>2 in section 3.1. The local charge neutrality 
of the TBBM is clearly inadequate to describe an ionic crystal for which a dominant part 
of the total energy is the Madelung sum of electrostatic pair terms. 10 A way to deal with 
this in tight binding was proposed by Majewski and Vogl 42 ' 59 based on a Hubbard-like 
hamiltonian of Kittler and Falicov. 60 In this scheme the total charge transfer at each site, 
A<7r, from (11) and (19) are taken as point charges. The hamiltonian is again 

H = H° + H' 

as in (22). Two terms make up H', the Madelung energy of the lattice of point charges and 
a positive energy that is quadratic in A<jr, namely [/rA<^; employing the well-known 
"Hubbard U" that acts to resist the accumulation of charge. This problem is solved self 
consistently. An extension of this scheme to allow the charge to be expressed as multipoles, 
not just monopoles, was proposed independently by Schelling et al. 61 and Finnis et al. 39. 
In the latter paper, the connection was made to density functional theory and the TBBM, 
so we will pursue the same argument here. As noticed by Elstner et al. 62 the Hohenberg- 
Kohn total energy in DFT can be expanded about some reference electron density, p°(r). 
If H° is the hamiltonian with effective potential generated by the reference density, and 
just as in section 2.1 its eigenfunctions are \n) then the total energy correct to second order 
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is 63 (e is the electron charge) 

£ (2) =£/n(n|if°|n) 



- j p (r)V^r)dr-E a H + E xc + E zz 



Mr) , Sb »(r')[. 



52J 

•6p(T)5p(r'y 

£"2 is the Hartree energy and £?° c and the exchange-correlation energy and potential 
belonging to the reference density, p°(r). The first two lines make up the Harris-Foulkes 
first order functional; we recognise the first line as the band energy, in fact the sum of 
occupied eigenvalues of the non self consistent input hamiltonian, and the second as the 
interaction term (double counting) plus the ion-ion pair potential, E zz . In the self consis- 
tent polarisable ion tight binding model 3S (SCTB) we approximate the last two lines by a 
generalised Madelung energy and a Hubbard energy, which adds a second order energy 5 
to (35) 

E 2 = \ e 2 Qn'L> B R 'L> rl Qrl + U ^R- ( 37 > 

These two terms represent the electron-electron interactions. All the exchange and cor- 
relation complexities are rolled into a single parameter, the Hubbard U. The first term 
in (37) is a classical interaction energy between point multipoles. The monopole term 
is just a straight forward sum of Coulomb energies, \e 2 A^rA^r/ / |R — R'|, while the 
generalised Madelung matrix is just the LMTO bare structure constant matrix (28), or to be 
precise -Br'l' rl = -(l/27r)(2^ + l)(2f + 1)-Br'l'rl- In general Qrl is the multipole 
moment of angular momentum L at site R. If we knew the charge density, which we don't 
in tight binding, then we could define the moment 



Qrl = j drp(r) r e Y L (r) 
for I > 0; while for £ = we'll have 



(38) 



Qr = Aq R r = \ ~r a, ?r- 

V 47T 

Although we don't know the charge density in tight binding, we know the eigenvectors of 
the hamiltonian and we can construct multipole moments from these. The monopole is of 
course proportional to the Mulliken charge transfer. Although in tight binding we don't 
even specify what the basis functions (3) are, we can take it that they comprise a radial part 
times an angular, spherical harmonic part, that is 

(r|RL) = fnt (|r - R|) Y L (r - R). (39) 

Then in terms of the eigenvector expansion coefficients (4), for I > we may define 



Qrl= ^ ^/ nC -" Ri ,4i»(Ri' 

L'L" n 



Qrl 



RL") (40) 
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in which the multipole moment operator is 

Qrl = r e Yl (r) , (41) 

which follows as a consequence of (38). If we expand out the matrix element of Qrl 
using (39) and (41) we have 

= Ai> £n e Cl'l"l, 

which introduces new tight binding parameters, A^£»^. Selection rules which are policed 
by the Gaunt coefficients (29) demand that there are only seven new parameters, or two if 
one has a basis of only s and p orbitals. These parameters are 



Aon 


= Aioi 


— A S pp 




Ana 


= Appd 


Ao22 


= A 2 02 


= Asdd 


A m 


= A 2n 


— Apdp 




A222 


= Addd 


A123 


= A213 


= A pd/ 




A224 


= Addg 



In fact these parameters are not entirely new, but are recognisable as the elements of crystal 
field theory — in the case I' = I" they are the quantities ( r e ). 65 ' 66 So it's perhaps not 
surprising that these new parameters introduce crystal field terms into the hamiltonian. 
These are off-diagonal, on-site terms that we have up to now taken to be zero. However 
they are crucial in describing the bands of, for example, the transition metal oxides as in 
figure 3. The generalised Madelung energy in (37) implies that the electrons are seeing an 
electrostatic potential due to the multipole moments at all the atomic sites. Indeed, if the 
electrostatic potential in the neighbourhood of the atom at site R is expanded into spherical 
waves, we could write, 

Vn(r) = J2 V *Lr e Y L (r) (42) 

L 

and using standard electrostatics the RL coefficient in this expansion is 

Vrl = ^ -Brxr/L' Qr'L'. 
R'L' 

Now in the same way that we arrived at (40), using (42) we can find the matrix elements 
of H' , namely 

#rl'rl" = A<7r. 5l'L" + e 2 ^2 Vrl c l<l»l- (43) 

L 

Now all the ingredients of the self consistent tight binding scheme are assembled. H° is 
given by its matrix elements, determined as in non self consistent tight binding, described 
in section 3. After solving the orthogonal, or non orthogonal eigenproblem and finding 



RL' 



Qrl 
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the eigenvector expansion coefficients, you build the multipole moments and using struc- 
ture constants find the components, Vrx, of the potential. Having also chosen the A and 
Hubbard U parameters, elements of H 1 are assembled and the eigenproblem is solved for 
H° + H' . This continues until self consistency. 
One or two extensions have been omitted here. 

1 . Only on-site matrix elements of H' are non zero in this self consistent scheme. In fact 
in the case of a non orthogonal basis, due to the explicit appearance of bond charge 
(see equation (11) and subsequent remarks) also intersite matrix elments of H' are 
introduced. This is important because it allows the hopping integrals themselves to be 
affected by the redistribution of charge, as might be intuitively expected. 6 ' 67 Details 
are to be found elsewhere. 5 23 

2. This scheme can be extended to admit spin polarisation in imitation of the local spin 
density approximation. This magnetic tight binding (figure 1) has also been described 
elsewhere and is omitted from these notes for brevity. 23 

Finally we should remark that the interatomic force is easily obtained in self consistent 
tight binding. Only the first and third terms in the TBBM (21) survive; in particular one 
still requires the derivatives of the matrix elements of H°. The only additional contribution 
to the force comes from the first term in (37); there is no contribution from the second term 
(or from the Stoner term in magnetic tight binding 68 ) because of the variational principle. 
Hence one requires only the classical electrostatic force on atom R, 

L 

which is consistent with the force theorem, 3134 and repairs the inconsistency of the band 
model mentioned in section 2.3. 

We illustrated the self consistent polarisable ion tight binding model (SCTB) in the 
study of phase transitions in Zr02 in section 3.1. It turns out that the extension of the 
point charge model to include polarisability introduces new physics that is essential in 
describing these phenomena. In particular the dipole polarisation of the anions drives the 
cubic to tetragonal transition. Furthermore, as seen in figure 3 the crystal field splitting 
of the cation d-bands is achieved naturally and the correct ordering is reproduced in cubic 
and octahedral crystal fields. Crystal field splitting is also largely responsible for the ligand 
bandwidth in the low symmetry rutile structure. 

4.1 Application to small molecules 

Now we will turn to a second example, the application to small molecules. The self consis- 
tent point charge model in this context and in the study of biological molecules has enjoyed 
enormous success thanks in particular to the work of Frauenheim, Elstner and colleagues. 69 
Here we demonstrate the SCTB model applied to the question of the polarisability 
of two small molecules, azulene and para-nitroaniline (pNA). Hopping parameters were 
taken from Horsfield et al. 70 and Hubbard U and A parameters chosen to to reproduce the 
ground state dipole moments predicted by the local density approximation. For azulene it is 
found that the self consistent point charge model is sufficient, but pNA cannot be described 
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Figure 6. Dipole moment as a function of applied electric field calculated using LSDA, solid lines, and SCTB, 
dotted lines. 71 LSDA calculations were made using a molecule LMTO program. 72 ' 73 The left hand figure shows 
the molecule azulene and the upper set of lines refer to the ground state and lower set to the so called Si ex- 
cited state. The right hand figure shows p-nitroaniline; the lower set are the ground state and the upper set the 
"zwitterionic" first excited state 



properly without dipole polarisability. 71 Figure 6 shows that the SCTB model provides a 
very accurate rendering of the dipole response to an applied electric field compared to 
LSDA calculations. We discuss now the two molecules in turn. 

1 . Azulene is a very interesting molecule having the same chemical formula as naptha- 
lene but comprising a five and seven membered ring instead of two six membered 
rings. According to Hiickel's "4n + 2 rule," a ring molecule is especially stable if it 
has N 7r-electrons and N = An + 2, where n is an integer. This is because this leads 
to a closed shell of 7r-electrons. 74 Hence benzene is stable, having n = 1. By a sim- 
ilar argument a seven membered ring has an unpaired electron which can be used to 
occupy an unpaired hole in a five membered ring. Hence the ground state of azulene 
possesses a large dipole moment. An excited state is created if the electron is returned 
to the seven membered ring. As shown to the left of figure 6 the ground state dipole 
moment is positive (the positive axis pointing to the right) while its sign is reversed in 
the first excited state. Here we use a device which is not quite legitimate, namely in 
both LSDA and SCTB an electron-hole pair is created and self consistency arrived at 
under this constraint. While a very crude approximation to an excited state 75 (given 
that LSDA is a ground state theory) this does provide a useful test of the validity of the 
SCTB model. Indeed it is quite remarkable how the SCTB faithfully reproduces the 
LSDA even to the extent of accurately reproducing the polarisability of both ground 
and excited states. (The polarisability is the linear response of the dipole moment to 
an applied electric field, namely the slope in these figures.) 

2. pNA is the archetypal "push-pull" chromophore. 76 In the ground state the dipole 
moment is small, but the first excited state is thought to be "zwitterionic," meaning 
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time (fs) time (fs) 

Figure 7. Charge transfer and bond current as a function of time in the relaxation of the S2 excited state in 
azulene. The upper panels show the excess charge on a "bridge" atom and on the rightmost atom in the seven 
membered ring (lower curve). The lower panels show the 7r-bond current in the "bridge" bond. 



that an electron transfers from the amine group on the right to the NO 2 group at the 
left increasing the dipole moment as shown on the right hand side of figure 6. Transfer 
of the electron through the 7r-system is called a push-pull process. Again the SCTB 
faithfully reproduces the LSDA with quantitative accuracy. We should mention again 
that it did not seem possible to obtain this result using a point charge self consistent 
tight binding model. 



4.2 Ring currents in azulene 

The SCTB model provides a simple scheme for the study of electron transfer as in the 
push-pull process. This is done by solving the time dependent Schrodinger equation using 
the hamiltonian H including electron-electron interactions. Indeed this is probably the 
simplest quantum mechanical model that goes beyond non interacting electrons. We have 
applied this approach to the relaxation of the S2 excited state in azulene with some quite 
spectacular results. 77 
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In terms of the density operator, the time dependent Schrodinger equation is 



±p = (ih)- 1 [H,p]-r(p-p ). 



We have added a damping term with time constant r^ 1 . This allows us to prepare the 
molecule in an excited state and relax it into the ground state whose density operator is po. 
The equation of motion is solved numerically using a simple leapfrog algorithm. While 
at the outset, the density matrix is real, during the dynamics it acquires complex matrix 
elements whose imaginary parts describe bond currents, 11. 



which is the total current flowing from atom R to atom R'. By selecting certain i-channels 
we can extract orbital contributions to j; in the present case of push-pull transfer we are 
interested in the current carried by the 7r-system of electrons. 

Figure 7 shows results of such a simulation in azulene, using a time constant r^ 1 = 
500 fs. Examine first the lower curve in the upper left panel. This is the excess total 
charge on the rightmost atom in the seven membered ring (see the inset in the top left 
panel). In the excited state, the dipole moment points to the left, that is, there is excess 
charge on this atom which transfers through the 7r-system to the left as the molecule relaxes 
into the ground state for which the dipole moment has opposite sign. The curve clearly 
show a smooth transfer of charge away from this site. However superimposed upon this 
is a series of oscillatory excursions in charge transfer, shown in a narrow time window 
by the lower curve in the upper right panel. Accompanying these oscillations are much 
larger fluctuations in the charge on the upper atom belonging to the "bridge" bond which 
is shared by both the five and seven membered rings. This excess charge is plotted in 
the upper curves of the upper left and right hand panels. As the upper and lower left 
hand panels show these oscillations die away, but analysis shows a quite characteristic 
frequency as seen in the right hand panels. The lower two panels show the 7r-bond current 
in the "bridge" bond. What is happening here is the setting up of ring currents in both 
rings whose directions are alternating with a period of a few femtoseconds. The ring 
currents at any one time are travelling in opposite senses in the two rings. This phenomena 
is a consequence of the electron-electron interaction, as we can verify by repeating the 
calculations using the non interacting hamiltonian, H°. Because two bonds enter each 
bridge atom but only one leaves, the opposing sense of the currents means that charge will 
accumulate on one of these atoms to the point at which the Coulomb repulsion (described 
by the Hubbard U) resists further current flow and indeed reverses its direction. Note that 
each current reversal (lower right panel) is mirrored by the alternating charge transfer on 
the bridge atoms (upper right panel). It is not yet understood what fixes the frequency at 
which the reversal happens or what it is that makes the molecule particularly susceptible 
to this instability. We note that these ring currents require a long lead-in time, on the order 
of the time constant, to become established and this is probably because the symmetry 
breaking comes about through numerical round-off in the computer. In a more detailed 
simulation coupling the electrons to the molecular vibrations, 78 this symmetry breaking 
will derive from the coupling. We can confirm that the great majority of the current is 
indeed carried by the 7r-electron system. 
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5 Last word 



The intention here has been to provide a practical introduction to the tight binding method 
and to motivate students to try it for themselves. While this is a long article it is mostly 
conspicuous for what is missing, rather than what is included. This is not surprising in 
view of the vast literature and considerable age of the tight binding approximation, but 
I've tried to bring out issues that are less widely discussed elsewhere. Regrettably no 
connection has been made to the semi empirical approaches in quantum chemistry that bear 
a close resemblence. This reflects the fact that physicists and chemists frequently discover 
the same science independently and often without much awareness of each other's work. 
I hope that some of the most glaring omissions will be covered by other authors in this 
volume. 8 ' 79 

Appendix 

Real spherical harmonics are described in ref [64]. One takes the conventional, complex 
spherical harmonics 80 and makes linear combinations to get the real and imaginary parts. 81 
Instead of m running from —I to I, m now runs from to I but for each m > 0, there 
are two real functions: Y£ which is (— l) m \/2 times the real part of Yi m ; and Y£ which 
is (— l) m \/2 times the imaginary part of Yg m . For m = 0, Yi m is anyway real, so we 
throw away Yf . We end up with the same number of functions, properly orthonormal. 
Specifically, 

Yf m = (-ir±(Y tm +Y tm ) 
Y/ m = (-ir^=(Y tm -Y lm ). 
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